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ABSTRACT 


Synthetic  seismograms  have  become  an  essential  tool 
in  the  interpretation  of  real  seismic  records.  Due  to 
the  low  cost  of  numerical  computation  and  the  possibility 
to  identify  and  study  individual  events  on  seismograms, 
asymptotic  ray  theory  (ART)  is  well  suited  for  this 
purpose. 

Asymptotic  ray  theory,  however,  has  its  limitations. 

It  can  neither  be  used  for  the  investigation  of  diffracted 
waves  in  shadow  zones,  nor  in  the  neighborhood  of  singular 
points  such  as  caustics  and  critical  points. 

The  aim  of  this  thesis  is  to  extend  the  applicability 
of  ART  to  include  the  effects  of  caustics.  First,  the 
essential  features  of  ART  are  presented.  Then,  using  the 
wave  method,  a  formal  integral  solution  for  an  arbitrary 
ray  in  a  layered,  vertically  inhomogeneous  medium  is 
derived.  Evaluating  this  integral  solution  using  a  modi¬ 
fied  version  of  the  third-order  saddle  point  approximation, 
a  new  expression  is  derived  for  the  amplitude  in  the 
vicinity  of  a  caustic.  This  result  is  shown  to  be  more 
accurate  than  the  amplitude  expression  currently  in  common 
use.  Pulse  distortion  due  to  caustics  is  also  investigated 
and  a  new  expression  is  derived  allowing  easy  determination 
of  the  phase  shift.  Finally,  numerical  results  (synthetic 
seismograms,  ray  diagrams,  time-distance  and  amplitude- 
distance  curves)  based  on  the  extended  ART  are  presented. 
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CHAPTER  1 


ASYMPTOTIC  RAY  THEORY  FOR  VERTICALLY 
INHOMOGENEOUS  MEDIA 

This  chapter  presents  the  basic  results  of  ART 
applied  to  the  problem  of  elastic  wave  propagation. 
Details  of  the  theory  can  be  found  in  Karal  and  Keller 
(1958) ,  Hron  and  Kanasewich  (1971) ,  and  Cerveny  and 
Ravindra  (1971). 

The  linearized  equations  of  motion  governing  small 
amplitude  wave  propagation  in  an  inhomogeneous,  per¬ 
fectly  elastic  and  isotropic  medium  are 

p-^-y  =  (X+y )  V  (V  *u)  +  yV2u  +  VX(V-u) 

at 

i.i 

+  Vyx(Vxu)  +  2(VyV)u 


where  u  is  the  particle  displacement  vector,  t  is  the 
time,  p  is  the  density  of  the  medium  and  X  and  y  its 
Lame  parameters. 

The  ray  method  assumes  that  the  complete  solution 
to  1.1  and  given  boundary  conditions  can  be  decomposed 
into  an  infinite  number  of  individual  contributions 
each  of  which  is  independently  a  solution  of  the 
problem.  It  is  further  assumed  that  an  individual 
contribution  can  be  expanded  in  an  asymptotic  series. 


m 


2 


called  the  ray  series,  in  inverse  powers  of  the  frequency 


u 


=  S(u))  e 


-10)T 


l 

k=0 


(iuj ) 


-k+ 

u. 


Here  S(co)  is  the  spectrum  of  the  wave  pulse,  t  and  u^ , 
which  are  functions  only  of  the  spatial  coordinates, 
are  called  the  phase  function  and  amplitude  coefficients 
of  the  ray  series. 

The  moving  surfaces  of  constant  phase,  t  =  x(x,y,z), 
are  called  wave  fronts  and  the  orthogonal  trajectories  of 

the  surfaces  rays.  In  general, any  lit  point  in  the  medium 
is  uniquely  defined  by  the  position  vector  X  =  5(T,q^,q2) 
where  the  curvilinear  coordinates  q^  and  q2,  called  the 
ray  coordinates , characterize  a  ray  while  the  phase  t 
defines  the  position  of  a  point  along  the  ray.  By  a  ray 
tube  we  will  mean  the  family  of  rays  the  parameters  of 
which  are  within  the  limits (q^ ,q^+dq1) and(q2 ,  q2+dq2) • 

On  any  given  wavefront  the  cross-sectional  area  da 
of  the  ray  tube  is 


da 


3Xy8X 
9q^X  3q2 


1.2 


Due  to  the  inverse  frequency  dependence  of  the  ray 
series,  it  is  apparent  that  the  first  term  is  the  most 
significant.  In  this  thesis,  only  the  first  term  of  the 


' 
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ray  series  will  be  considered.  This  leads  to  the  zeroth- 
order  solution  or  zeroth  approximation  of  ART: 


u  =  S  (to) 


-iwT 

e 


U (x,y , z) 


Furthermore,  it  suffices  to  consider  only  an  arbitrary 

ray  as  similar  methods  can  be  applied  to  compute  all 

other  rays  composing  the  total  solution. 

In  vertically  inhomogeneous  media  all  rays  are  plane 

curves.  Thus  a  Cartesian  coordinate  can  be  adopted  so 

that  the  entire  ray  is  contained  in  the  vertical  xz-plane. 

Then  three  orthogonal  unit  vectors  n  ,  nc__  and  nnT1  can  be 

p  bv  SH 

introduced  at  every  point  of  the  ray  such  that  n^  is 
tangent  to  the  ray  pointing  in  the  direction  of  propaga- 


tion;  ngv  is  perpendicular  to  the  ray  and  lies  in  the 
vertical  plane  such  that  its  projection  on  the  x-axis 
is  positive;  and  ngH  is  a  unit  vector  in  the  direction 
of  the  y-axis. 


Substituting  the  zeroth-order  solution  into  the 
equations  of  motion  one  finds  that 


u 


c  ,  v  “1C0T  TT 

S  (to)  e  U .  n  . 

1  1 


where  ni  may  be  n  ,  ngv  or  ngH,  and  are  the 

corresponding  amplitude  coefficients.  This  means  that 


'  r 


- 


in  the  zeroth  order  approximation  of  ART  three  kinds  of 
elastic  waves,  P  waves,  SV  and  SH  waves,  can  propagate 
within  vertically  inhomogeneous  media. 


In  all  three  cases,  the  phase  function  t  satisfies 
the  well  known  eikonal  equation  (see  e.g.  Cerveny  and 
Ravindra  1971,  eq.  2.15  and  2.16)  which  can  be  solved  to 
yield 


t  (M) 


t(Mq)  + 


ds 

v(x,y , z) 


1.3 


where  s  is  distance  measured  along  the  ray;  M  and  Mq  are 
points  on  the  ray;  v  is  the  local  wave  velocity  given  by 


for  P  or  compressional  waves  and 


for  SH  and  SV  or  shear  waves.  The  ray,  everywhere 
parallel  to  the  gradient  of  the  phase  function  t,  is 
obtained  by  solving  the  Euler  equations  for  the  extremals 
of  the  integral  in  1.3. 

The  amplitude  coefficients  are  governed  by  the 
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transport  equations  (see  e.g.  Cerveny  and  Ravindra  1971 
eq.  2.25 1  and  2.37)  which  can  be  integrated  to  give 


Ui(M)  =  (-1) 


m 


(vpdcr) 


M 


(vpdcr) 


M 


U. (M  ) ;  i=l ,2,3 
1  o 


1.4 


L  -J 

Here  the  index  m  is  zero  for  all  P  and  SH  rays.  For  SV 
rays  it  assumes  the  value  one  if  the  ray  has  a  turning 
point  between  Mq  and  M.  Otherwise  m  equals  zero. 

The  amplitude  coefficients  1.4  have  a  very  simple 
physical  interpretation.  The  magnitude  of  UL  varies 
precisely  in  accordance  with  the  principle  of  conserva¬ 
tion  of  energy  flux  within  a  narrow  tube  of  rays  where 
the  energy  propagates  only  along  the  rays  and  no  energy 
flows  through  the  side  walls  of  the  ray  tube.  Thus  each 
ray  is  a  distinctive  path  of  energy  transport  between 
the  source  and  receiver. 

The  considerations  so  far  require  that  the  elastic 
parameters  of  the  medium  be  continuous  along  the  ray. 

The  results  can  be  easily  generalized  to  include 
reflections  and  transmissions  at  interfaces.  It  is 
assumed  that  the  medium  consists  of  layers  in  welded 
contact  so  that  displacement  and  traction  are  continuous 
across  the  interfaces.  From  the  six  equations  express¬ 
ing  this  continuity,  it  follows  that  in  the  zeroth 
approximation  of  ART  (see  e.g.  Cerveny  and  Ravindra  1971, 
section  2.3) 
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1.  SH  waves  propagate  in  complete  independence  of  the 
other  waves. 

2.  conversion  between  P  and  SV  waves  can  occur  by 
reflection  and  transmission  at  interfaces. 

3.  energy  partitioning  at  interfaces  can  be  described 

by  plane  wave  reflection  and  transmission  coefficients. 
Using  1.4  for  the  continuation  of  the  displacement 
vector  along  the  ray  within  a  continuous  layer  and  using 
plane  wave  reflection  and  transmission  coefficients  for 
continuation  across  interfaces,  it  is  easy  to  construct 
the  displacement  vector  for  a  ray  that  has  encountered 
an  arbitrary  number  of  interfaces  between  the  points  Mq 
and  M  (see  Fig.  1) .  If  we  write 


u(M)  =  S(w)U(M)  e”la)Tn(M) 


1.5 


where  n(M)  may  be  n  (M)  ,  n_,TT(M)  ,  or  n^CM) 

p  o  V  on 

the  type  of  wave  observed  at  M,  then 


depending  on 


k+1 

t (m)  =  t(m  )  +  y 
°  j=i 


o 


°j-i 


ds 

v 


?  0  EM,  0,,.  EM  1.6 
o  o  k+1 


and 


p 

(VP)M 

* 

(vp)o7 

(-1) 

o 

k 

3 

L  (M) 

(VP)M 

1  1 

(vp)  + 

j=l 

3 

R(0  .) 


U  (M) 


1.7 


- 

■ 
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Figure  1 


Typical  ray  path  in  a  layered,  vertically 
inhomogeneous  medium.  and  M  are 

respectively  the  source  and  receiver 
points;  CK  ,  j=l,...,k  are  the  points 
of  incidence  at  interfaces. 
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X 


where,  for  simplicity,  U(Mq)  has  been  set  to  unity  (only 
relative  amplitudes  are  significant  in  synthetic  seismo¬ 
gram  computations) ;  e  is  the  number  of  turning  SV  ray 
segments;  R(CK)  is  the  appropriate  reflection  and  trans¬ 
mission  coefficient  at  the  point  of  incidence  0  ^  ;  the 
superscripts  -  and  +  indicate  properties  at  the  point  CK 
on  the  side  of  the  incident  wave  and  on  the  side  of  the 
reflected  or  transmitted  waves,  respectively;  and  the 
function  L (M) ,  called  the  spreading  function,  is  given  by 


L  (M)  = 


da  (M) 

^  k 

da (ot) 

D 

da  (Mq) 

TT 

3=1 

da  (0  . ) 

3 

h 


In  general,  the  spreading  function  may  be  evaluated 

using  1.2.  However,  for  a  point  source  in  a  vertically 

inhomogeneous  medium,  a  rather  simple  expression  may  be 

obtained.  We  assume  that  in  the  immediate  vicinity  of 

the  source  the  medium  is  essentially  homogeneous.  Let 

M  be  a  point  at  unit  distance  from  the  source  within 
o 

this  homogeneous  region.  Then  since  the  wavefront  is 


spherical  at  Mq : 


da  (M  )  =  sine  d0  d<J> 

G  O  O  O 

where  6q  and  4>q  are  spherical  polar  coordinates  des¬ 
cribing  the  ray  at  M  .  From  simple  geometrical  considera¬ 
tions  the  ray  tube  cross-sectional  area  at  M  can  be 


s 
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written  as  (see  Fig.  2) 


da  (M)  =  r-l~  cos9(M)d(}>  d0  . 

do  O  O 


1.8 


o 


Furthermore,  at  the  point  CK  at  the  end  of  the  j-th  ray 
segment  we  have  (see  Fig.  2) 


cosQ  (0+) 


da  (0  . ) 
D 


cosO (0 j ) 


so  that  the  spreading  function  can  be  written  as 


3r cosO (M) 

% 

k 

cos0 (0+) 

1 

r30  sin0 

L  o  o 

n 

j=] 

. 

cosO (0 j ) 

Equations  1.5  to  1.7  and  1.9  constitute  the 
zeroth  approximation  of  the  particle  displacement 
provided  by  ART  in  its  basic  form.  Its  deficiency  is 
apparent  as  an  infinite  amplitude  is  predicted  at  points 
where  the  ray  tube  cross-sectional  area  da  vanishes. 
Since  da  is  also  the  Jacobian  of  transformation  between 
spatial  and  ray  coordinates,  it  is  clear  that  the  above 
results  are  invalid  not  only  at  these  points  but  also 
in  regions  along  the  ray  beyond  these  points  (Babich 
19  61)  .  The  following  chapters  aim  to  resolve  these 
difficulties . 
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Figure  2  Ray  tube  cross-sectional  area. 


(a)  Cross-section  at  the  ooint  M  of  an 
elementary  ray  tube  originating  at  the 
point  M. 

(b)  Change  of  the  cross-sectional  area 
across  an  interface. 
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(D) 
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CHAPTER  2 

RAY  AMPLITUDE  BY  WAVE  METHOD 

A.  Introduction 

In  the  previous  chapter  the  propagation  of  small 
amplitude  elastic  waves  was  studied  using  ART.  The 
theory  is  useful  but  unfortunately,  under  certain 
circumstances  the  results  are  invalid  and  more  exact 
methods  must  be  used. 

In  this  chapter  a  formal  integral  solution  to  the 
problem  of  a  point  source  in  a  plane-layered,  vertically 
inhomogeneous  elastic  medium  will  be  obtained.  In  the 
following  chapters,  this  integral  solution  will  be 
evaluated  to  provide  an  extension  to  ART. 

B.  The  Elastodynamic  Equations 

We  consider  a  perfectly  elastic,  isotropic  medium 
in  which  the  wave  velocities  and  density  are  functions 
only  of  the  vertical  coordinate  z.  For  symmetry  in  the 
solution,  we  consider  an  explosive  point  source  located 
at  the  point  with  cylindrical  coordinate  (r,z,(J>)  = 
(0,zq,0)  and  described  by  the  source  function 

6 (z-z  ) 6 (r) 

P  (r,  z,  t)  =  PQ  (t) 


2i\r 


Since  the  ray  path  in  a  vertically  inhomogeneous  medium 
is  a  plane  curve,  the  azimuthal  angle  cf>  may  be  ignored. 
Within  the  medium,  the  equations  of  motion  are 
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f 


2.1 


where  the  body  force  per  unit  volume  f  is  derivable 
from  the  source  function  by 


and  the  stress  dyadic  o  is  related  to  the  particle 
displacement  u  via  the  constitutive  equations 


cr  =  (XV-u)I  +  y  Vu+  (uV ) T 


I  is  the  unit  dyadic  and  the  superscript  T  indicates  a 
transpose . 

The  partial  differential  equations  2.1  may  be 
considerably  simplified  if  transformed  coordinates  are 
used.  We  will  use  the  Fourier  transform  in  time  defined 
by 


oo 


2.2 


—  00 


. 
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and  the  zeroth-order  Hankel  transform  in  the  cylindrical 
radius  r  defined  by 


U  (  p  ,  Z  ,  (jl)  ) 


u(r , z,co) 


OQ 

u  (r ,  z ,  to)  r  Jq  (topt )  dr 
o 


2 

co 


f  CO 

u(p,z,oj)  pJQ  (wpr)  dp 
o 


For  simplicity  of  notation,  a  function  is  distinguished 
from  its  transform  only  by  its  arguments. 

With  these  transformations,  the  equations  of  motion 
and  constitutive  equations  are  reduced  to  a  system  of 
first  order  ordinary  differential  equations 


dV  . 

3—  =  loaAV  +  w 
dz 


where 


(p,  Z,0)) 


zi  !_ 

pr  9r 


(rur) , 


ium  , 
z 


-1  9_ 

iwpr  9r 


2.4 


are  the  doubly  transformed  components  of  displacement  and 
stress ; 


w 


(p,t,w) 


_1 _  3_ 

iwpr  8r 


1 


(rf  )  1 
r 

J 


2.5 


are  the  doubly  transformed  force  functions;  and  the 
matrix  A  is  given  by 
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O 


0 


P 


1 

y 


o 


o 


p 


p 


A 


£X 


1 


O 


o 


X+2]i 


A+2 


P 


0 


o 


X+2\i  A.+2y 


The  differential  system  2.3  is  well  known  in  seismo- 
logical  literature  and  arises  also  in  the  study  of  SH  and 
acoustical  waves.  After  Gilbert  and  Backus  (1966),  we 
will  call  the  non-singular  matrix  Vf  a  fundamental  matrix 
of  this  differential  system  if  it  is  a  solution  of  the 
homogeneous  form  of  the  system,  i.e.  if  each  of  its 
columns  satisfies 


2.6 


If  a  fundamental  matrix  has  the  property  that  Vf(zQ)  the 
identity  matrix  then  it  is  called  a  propagator  from  zq 
and  is  denoted  by  P(z;z  ).  Propagators  can  be  constructed 
from  any  fundamental  matrix  by 


=-l 

P  ( z ;  z  )  =  V,-(z)  V.  (z  ) 
o  t  to 


2.7 


, 


» 


17 


Using  propagators,  the  solution  of  the  homogeneous 
system  2.6  can  be  easily  expressed  as 


V(p, z,w)  =  P (z 


zq)  V(p,zQ,(jo) 


2.8 


while  that  for  the  inhomogeneous  system  2.3  is 


r  z 


V(p,z,a>)  =  P(z,-zq)  V(p,z  ,a>)  + 


P(z;£)w(p,£,u>)d£  .  2.9 


Thus  we  see  that  the  basic  problem  in  solving  the 
elastodynamic  equations  is  to  find  a  fundamental 
matrix  V^. 

C.  The  Fundamental  Matrix 

The  matrix  A  has  eigenvalues ±q  and±q  where 

p  s 


v  and  v  being  the  local  P  and  S  wave  velocities.  The 
p  s 

corresponding  eigenvectors,  normalized  with  respect  to 
energy  flux  in  the  z-direction,  are 


V+  =  - i - jr  (-p,yr,  +q  ,  ±2ypq  ) 

±P  (2pq  ) ^  P  P 

ir 


v+  =  - - - -  (±q  ,  ±2npq  ,  -p,  -uD 

±S  (2pqs)^  S  S 


Sr 


. 
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where 


If  we  define  a  matrix  N  using  the  eigenvectors 
as  columns 


then,  using  N  ,  we  can  diagonalize  the  matrix  A  by  a 
similarity  transformation  yielding  the  matrix  of 
eigenvalues 


0 


0 


0 


Q 


A  N  = 


0 


0 


0  0  -qs  0 


0  0  0 


Thus  it  follows  that 
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N-1  A  =  q  n  1 

Using  this  relationship  the  homogeneous  system  2.6 
can  be  written  as 

(I"1  =  iooQCN-1  $)  +  BCN-1 

dz 


where 


B  = 


dN 


-1 


dz 


N 


Notice  that  B  involves  the  derivatives  of  the  elastic 

parameters  as  well  as  terms  in  inverse  powers  of  q  and 

q  ;  also  it  is  one  order  in  frequence  lower  than  iwQ  . 
s 

Therefore,  if  the  problem  is  such  that  the  WKBJ  conditions 

are  satisfied,  i.e.,  if  the  variations  in  the  elastic 

parameters  of  the  medium  are  small  over  a  wavelength  and 

that  the  solution  is  not  required  near  turning  points 

where  either  q  or  q  vanishes,  then  we  can  ignore  the 

P  s 

term  in  B  and  approximate  the  homogeneous  system  2.6  by 


d_ 

dz 


(N  1  V) 


iooQ  (N 


id) 


z 


QU)<H 


for  which  e 


—  00 


is  a  solution  matrix. 


H 


' 


'  • 
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In  the  WKBJ  approximation,  therefore,  a  fundamental 
matrix  of  the  elastodynamic  equations  2.3  is 


iw 


QU)d£ 


Vf  =  N  e 


2.10 


D.  The  Ray  Expansion 

Cisternas  et  al  (1973)  and  Kennett  (1974)  have 
shown  that  for  a  stack  of  plane-parallel,  homogeneous 
layers  the  solution  to  the  elastodynamic  equations  can 
be  decomposed  into  an  infinite  number  of  individual 
contributions.  Identical  with  the  rays  in  ART,  each  of 
these  contributions  is  independently  a  solution  of  2.3 
and  characterized  by  a  distinctive  path  of  energy 
transport.  For  media  consisting  of  vertically  inhomo¬ 
geneous  layers  where  the  WKBJ  approximation  is  valid, 
the  ray  expansion  is  again  valid  (Kennett  1974) . 

In  the  following,  the  complete  wave  solution  to 
2.3  for  a  layered,  inhomogeneous  medium  need  not  be 
attempted.  For  comparison  with  ART  results,  it  suffices 
to  consider  the  contribution  of  an  arbitrary  ray  in 
the  ray  expansion. 

E.  The  Source  Wave 

We  will  denote  as  the  source  region  a  volume  about 
the  source  the  linear  dimension  of  which  is  small  enough 
so  that  the  elastic  parameters  within  it  can  be  regarded 
as  constant  but  large  enough  compared  with  the  wavelength. 


. 


' 


\ 
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For  weakly  inhomogeneous  media  where  the  WKBJ  conditions 

are  satisfied  such  a  region  can  always  be  defined. 

By  a  source  wave  we  will  mean  a  wave  propagating 

within  the  source  region  from  the  source  without 

experiencing  any  reflection  or  transmission  at  interfaces. 

We  assume  that  both  the  source  and  receiver  are 

located  on  a  free  surface  at  z  =  z  so  that  waves  exist 

o 

only  for  z<zq.  Thus  the  problem  for  a  source  wave  is 
equivalent  to  that  of  wave  propagation  within  a  homo¬ 
geneous  half  space  for  which  2.9  provides  the  complete 
solution.  Using  2.5,  2.7,  2.10  and  the  radiation 

condition  this  solution  is 


VT(p,Ms,oj) 


icoP  (tjo)e 
o 


—  1GJ 


r  z 


z  V5)ds 

o 


"p*2,,!  (2py  ms(2pVm0 


(p,  -pr,  -qp,  2ypqp) ^ ;  zQ>zs 


2.11 


where  M  = ( r  ,z  )  denotes  a  point  within  the  source 
s  s '  s 

region. 

For  the  computation  of  synthetic  seismograms  the 
particle  displacement  u  is  desired.  This  can  be  obtained 
from  the  verticle  displacement  u  by 


u 


u 

z 


(n-nz) 


n 


where  n  is  the  unit  vector  in  the  positive  z-direction 
z 


and  n  may  be  n  or  n-,..  for  P  and  SV  waves  respectively . 
-1  p  bv 

-> 

Using  this  relation  and  the  definition  2.4  of  V, 
the  source  wave  displacement  at  M  is 

b 


-*•  t  »»  \  S  ((a) ) 

u(p,M  ,0))  =  - - y 

s'  iO)q  (M  ) 

p  o 


3* 

Cp/q  )  2cos0 


M 


o  -10) 
e 


qpU>d5 


n 


(p/qD) ^cose 


M 


2.12 


where  0  is  the  acute  angle  between  n  and  n  and 

p  z 


S(03)  =  ! 


ia)po  (oi)  qp 


I  2  TT  P  Vp  COS0 


i  M 


_J 


is,  apart  from  a  multiplicative  constant,  the  spectrum  o 
the  source  pulse. 


F.  Propagation  Along  a  Ray 

Outside  the  source  region  wave  propagation  is  des¬ 
cribed  by  the  homogeneous  system  2.6.  Independent 
solutions  are  columns  of  the  fundamental  matrix  : 


io) 

e 


Z  ±q_U)d5 

CO  i - 


representing  upgoing  and  downgoing  P  waves  and 


io) 

e 


z 


±q 


—  oo 


S 


U)d£ 


representing  upgoing  and  downgoing  SV  waves. 


C ,  and 
±P 


i 


^  o 


C+  are  arbitrary  constants. 

—  s 

At  any  two  points  M^=(r^,z^)  and  M2=(r 2,Z2^  within 
a  continuous  volume  outside  the  source  region,  the  wave 
solutions  are  related  through  the  propagator  by  2.8,  i.e. 


QC5)d£  , 

N  (M1)V(p,M1,tu)  2.13 


V  (p  ,M^  •  co)  =  N(M2)  e  z1 


Assuming  turning  points  do  not  exist  along  the  ray 
so  that  the  WKBJ  approximation  is  valid  throughout  the 
region  under  consideration,  2.13  together  with  2.11 
can  be  used  to  show  that  if  the  wave  at  is  of  a 
given  type,  a  downgoing  P  wave  say,  then  it  will  also 
be  of  the  same  type  at  M2 .  In  other  words,  reflections 
off  the  velocity  gradient,  and  their  associated  mode 
conversions,  are  negligible  in  the  WKBJ  approximation. 
The  complex  amplitudes  at  the  two  points  are  related  by 


2.14 


where,  for  P  waves 


q  =  q 


P 


2.15 


and  for  SV  waves 


2.16 


I  m'  l 
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and 


u(p,M,w)  =  u(p,M,oo)n 


If  a  turning  point  exists  along  the  ray  between  M^ 
and  M2 ,  the  WKBJ  result  2.14  is  no  longer  valid.  However, 
only  a  slight  modification  is  needed  to  provide  the 
correct  expression.  Appendix  1  shows  that  the  effect  of 


a  turning  point  can  be  represented  by  the  reflection 

.  TT 

coefficient  e  •  Suppose  a  turning  point  exists  at  (r*,z*), 
then  the  complex  amplitudes  at  and  M2  are  related  by 


TT 


f  f Z 


G(M,  )  X2  "  la3\ 

r  „  »  ,  ,  »  m  1  C 

u (p ,M  , oj )  =  (-1)  G^M  ^  e 


q(5)  d^ 


+ 


r  z. 


2.17 


q(Od£ 


•  u (p,M1 ,w) . 


The  (-l)m  factor  is  introduced  for  the  same  reason  as  in 
Chapter  1  for  SV  waves. 


G.  WKBJ  Solution  for  an  Arbitrary  Ray 

When  a  wave  is  incident  upon  an  interface,  it  will 
either  be  reflected  or  transmitted  and  mode  conversion 
can  occur.  The  complex  amplitude  u^.  of  the  incident 
wave  and  the  complex  amplitude  uD  of  the  reflected  or 
transmitted  wave  are  related  through  the  appropriate 
plane  wave  reflection  and  transmission  coefficients  R 
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(as  given  by  Cerveny  and  Ravindra  1971,  section  2.3): 


2.18 


Ru 


u 


I 


R 


Thus,  starting  with  the  source  wave  displacement 
2.12,  the  relations  2.14,  2.17  and  2.18  can  be  used  to 
construct  the  particle  displacement  of  an  arbitrary  ray 
which  has  gone  through  n  turning  points  and  suffered  k 
reflections  and  transmissions  at  the  points  0 y 
j  =  1 , 2 , . . . ,k ,  before  reaching  receiver  point  M.  The 
particle  displacement  at  M  is 


S  (co )  (-l)£G(Mo) 


-iojip  (p)  tin?-  ,  G(0.) 


3—  R  (0  .) 
3 


u(p,M,<d)  ia)  )G(M) 

p  o 


2.19 


where 


k+1 

*Kp)  =  l'  ^  a  (P) 
j=l  : 


with 


z  . 


.  (p)  =  J  q  (p,  V  d^ 

3  ir? 


if  the  ray  has  no  turning  point  between  0 j j _i ' z j -i 
and  Oj=(rj,Zj),  and 


rz 


j 


ipj  (p)  =  2 


q(p  rO  d? 


* 
z  . 
3 


{■ 


. 


- 


if  the  ray  has  a  turning  point  at  (r^,z*)  between  0j_^ 

and  0  .  . 

D 


Inverting  the  Hankel  transform  and  keeping  only  the 
first  term  of  the  asymptotic  expansion  for  the  Hankel 
function  (Abramowitz  and  Stegun  196  5)  ,  we  obtain  from 
2.19  the  solution  in  the  frequency  domain 


u(M,oo)  =  (-1)  S  ( 03 ) 


0) 


27rr 


\  i  (n5-  -  t) 


ftp)  e-i“^P'r>dp 


2.20 


where 


f  (P)  = 


P^G(Mq) 
q  (Mo)G(M) 


k 

rr 

j=i 


0(0^) 

G(0j) 


R  (0  .  ) 
j 


and 


2.21 


x(p,r)  =  pr  +  iKp)  2.22 

is  the  phase  function.  This  is  the  plane  wave  approxi¬ 
mation  of  the  WKBJ  solution  for  the  displacement  amplitude 


of  an  arbitrary  ray. 


. 
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CHAPTER  3 

PHASE  SHIFT  DUE  TO  CAUSTICS 
A.  Introduction 

A  caustic  may  be  defined  as  a  point  at  which,  dcr , 
the  cross-sectional  area  of  the  ray  tube  vanishes.  The 
basic  theoretical  problems  connected  with  caustics  have 
been  investigated  by  Brekhovskih  (1960) ,  Ludwig  (1966) , 
Sato  (1969)  and  others. 

Essentially,  in  the  high  frequency  approximation, 
the  dynamic  properties  of  the  wave  field  are  character¬ 
ized  by  a  tt/2  phase  shift  of  the  harmonic  components  at 
the  caustic  and  a  concentration  of  wave  energy  in  its 
immediate  neighborhood.  For  pulse  waveforms,  this  tt/2 
phase  shift  results  in  a  marked  distortion  of  the  pulse 
shape,  providing  a  useful  criterion  for  proper  identifica¬ 
tion  of  later  phases  on  a  seismogram  (Choy  and  Richards 
1975) .  Similarly,  the  larger  amplitude  near  a  caustic 
renders  it  important  for  the  interpretation  of  seismo¬ 
grams  . 

In  this  chapter,  the  phase  shift  at  a  caustic  will 
be  established.  This  will  lead  to  a  useful  expression 
relating  the  phase  shift  to  the  number  of  turning  points 
of  the  ray.  Ray  amplitudes  near  a  caustic  will  be  dealt 
with  in  the  next  chapter. 

B.  Relation  between  Turning  Points  and  Caustics 

The  occurrence  of  caustics  and  their  associated 


' 
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phase  shifts  is  closely  related  to  the  occurrence  of 
turning  points  of  totally  reflected  waves.  In  fact,  the 
pulse  distortion  observed  in  totally  reflected  waves  had 
often  been  erroneously  ascribed  to  turning  points.  The 
origin  of  this  mistaken  association  is  probably  due  to 
the  fact  that  turning  points  coincide  with  caustics  in 
the  limiting  case  of  totally  reflected  plane  waves 

(Tolstoy  1968,  Silbiger  1968). 

In  general,  for  waves  generated  by  a  point  source, 
caustics  are  distinct  from  turning  points.  Indeed,  a 
pulse  may  travel  through  a  turning  point  and  yet  not 
form  a  caustic  and  therefore  not  experience  any  distor¬ 
tion.  Direct  P  waves  in  the  Earth  are  obvious  examples 
of  this  situation. 

A  caustic  corresponds  to  the  vanishing  of  the  ray 
tube  cross-sectional  area.  For  vertically  inhomogeneous 
media,  equation  1.8  implies  that  3r/9p  vanishes  at  a 
caustic,  where  p  =  sin  0/v  is  the  geometrical  ray  para¬ 
meter.  Thus  the  existence  of  caustics  can  be  easily 
established  by  examining  the  variations  in  the  sign  of 
the  derivative  8r/3p  along  the  ray  path. 

For  this  purpose,  it  is  convenient  to  divide  the  total 

ray  path  into  segments  containing  single  turning  points 
(called  here  turning  segments) ,  and  segments  with  no 
turning  points  (called  here  non-turning  segments) .  We 
will  denote  by  c^  the  ray  segment  bounded  by  s^_^  and  s^, 
where  s  is  the  distance  along  the  ray  measured  from  the 


. 
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/  9rv 

source  point.  Then,  the  value  of  on  the  raY  at  s 

is  given  by  the  sum  of  contributions  from  individual 
segments : 


For  non-turning  segments,  its  contribution  is  given 
by  the  standard  integral 


(  9r\ 

3p  j 

c  . 

i 


max  ( z^_^  ,  z^) 


min ( z  .  , , z . ) 
l-l  l 


d£ 

2  3 
v  q 


>  0 


3.1 


where  z .  and  z .  are  the  z-coordinates  of  the  points 
i-l  l 

denoted  by  s.  .  and  s.  . 

i-l  l 

For  turning  segments,  more  care  is  required.  We 
assume  that  the  velocity  function  near  a  turning  point  can 
be  expanded  in  a  Taylor  series  in  powers  of  z- z*  ,  z* 
being  the  z-coordinate  of  the  turning  point.  We  can  always 
define  the  turning  segments  short  enough  that  only  the 
linear  term  in  the  expansion  need  be  considered  in  approxi¬ 
mating  the  velocity  function  over  the  entire  segment;  i.e. 


v(z)  =  v(z*)  -  g(z-z*) 


3.2 


where 


dv(z*) 

dz 


g 


Js£ 


is  assumed  to  be  non-zero.  This  assumption  is  invalid  only 

in  the  special  cases  where  the  turning  point  coincides  with 

an  extremum  of  the  velocity  function  and  higher  order  terms 

must  be  considered.  Using  3.2,  the  contribution  of  a  turn¬ 
er 

ing  segment  to  can  easily  be  found  to  be 


/  3r 
'  9P 


ci 


-1 

2 

gp 


\ 


/  l 


vq 

z  .  -i 
l-l 


l  vq 


<  0 


3.3 


We  now  consider  a  ray  with  n  turning  points  at  s?  , 

i=l ,2 , . . . ,n.  The  number  of  caustics  associated  with  this 

ray  can  be  found  by  simply  counting  the  number  of  sign 

changes  in  (x— )  as  s  varies  from  0  to  s.„ ,  sA.  being  the 

^  3p  s  MM 

distance  of  the  observation  point  along  the  ray. 

3 

Between  the  source  and  the  first  turning  point  (g--)  g  has 

contributions  only  from  non-turning  segments  and  therefore 

3  IT 

is  always  positive.  As  s  approaches  s* ,  s 

tends  to  positive  infinity  since  the  integrand  in  3.1  has  a 
non-integrable  singularity  at  the  turning  point  where  q 
vanishes . 


At  the  turning  point  (|^)s  is  undefined.  Physically, 
this  is  expected  since  a  ray  with  a  slightly  larger  ray 
parameter  than  a  reference  ray  can  not  possibly  penetrate 
to  the  same  depth  as  the  turning  point  of  the  reference  ray. 
In  this  case,  the  cross-sectional  area  of  the  ray  tube  can 
be  expressed  alternatively  by 


J 
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da  =  r 


-H-|d0  dcj) 

89  1  o  o 


o 


Since  8z/86  is  obviously  positive  at  a  turning  point,  da 


o 


is  finite  and  caustics  never  coincide  with  turning  points 
for  waves  from  point  sources.  Therefore  there  can  be  no 
caustics  before  the  first  turning  point. 

Just  past  the  turning  point,  for  s  =  s*  +  6  where  6 

3  IT 

is  positive  and  infinitesimally  small,  (-^) g  jumps 
discontinuously  to  negative  infinity  by  3.3.  In  general 


9_r  \ 
8p  i 

\  I 


->  ±  oo  as  s  -*  s' 


< 

> 


s* 


3.4 


Between  consecutive  turning  points  s|  and  s*+^ 

additional  contributions  again  come  only  from  non-turning 

3  T* 

segments.  Therefore  (^—)  increases  monotonically  from 

o  p  S 

negative  to  positive  infinity  vanishing  exactly  once  for 

some  s  where  s*  <  s  <  s*  , .  In  other  words,  exactly  one 
c  1C  1+1 

caustic  must  exist  between  two  consecutive  turning  points. 

Between  the  last  turning  point  s*  and  the  observation 

point  s^  an  additional  caustic  may  or  may  not  exist.  Using 

3  ^ 

3.4  and  keeping  in  mind  that  (-^ — )  increases  monotonically 

dp  S 

from  s*  +  6  to  s„,  we  conclude  that  an  additional  caustic 
n  M 

exists  if  and  only  if  >  0. 


8p  s 


M 


The  above  considerations  are  summarized  schematically 

3  IT 

in  Fig.  3  where  )  is  plotted  as  a  function  of  s.  It  is 

dp  S 

clear  that  for  a  ray  with  n  turning  points,  the  number  of 
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caustics  formed  along  the  ray  is  given  by 

nc  =  n  -  *5  [l-sgn(||)]M  3.5 

Ray  diagrams  illustrating  this  relationship  for  layered 
media  are  shown  in  Fig.  4. 

C.  Phase  Shift  at  a  Caustic 

It  is  well  known  that  the  high  frequency  harmonic 
components  of  a  wave  experience  a  discontinuous  tt/2  phase 
shift  as  the  wave  propagates  through  a  caustic  along  a 
ray.  This  tt/2  phase  shift  has  been  derived  by  Landau  and 
Lifshitz  (1951),  Kay  and  Keller  (1954),  among  others.  In 
this  section,  an  alternative  derivation  is  presented.  It 
not  only  produces  the  expected  phase  shift  but  also  pro¬ 
vides  a  rule  for  extending  the  validity  of  the  asymptotic 
ray  theory  solution  to  regions  beyond  a  caustic. 

We  assume  that  a  caustic  is  formed  at  the  point 

M  =  (r  ,z  )  .  Analogous  to  the  source  region  of  the  previous 
c  c  c 

chapter,  we  define  the  caustic  region  as  a  neighborhood  of 

M  small  enough  so  that  it  can  be  regarded  as  homogeneous 
c 

but  large  compared  with  the  wavelength.  To  study  the 
effects  of  passage  through  a  caustic,  it  suffices  to  con¬ 
sider  wave  propagation  within  the  homogeneous  caustic  region. 

Inside  the  caustic  region,  independent  solutions  of  the 
transformed  wave  equations  2.6  are  exp(±icoqz). 

Therefore  the  general  solution  in  the  frequency  domain  is 


Figure  3 


Variation  of  3r/3p  with  arc  length  s  along 
a  typical  ray  path.  Dependence  of  n  on  n 


and  (3r/3p)M  is  apparent. 
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Figure  4  Examples  of  the  application  of  equation  3.5. 

The  ray  with  the  smaller  (bigger)  ray  para¬ 
meter  is  labelled  ray  1(2).  Solid  (dashed) 
curves  represent  P  (SV)  ray  segments. 
Caustics  are  marked  by  X. 
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u 


(r,z,w)  =  l 


od  a  (p) 


-ico  (pr±qz) 


3.6 


(r) 


dp 


where  the  summation  is  over  the  ±  indices;  and  a+ (p)  are 

functions  to  be  determined.  We  prescribe  the  initial 

conditions  at  a  point  M  =  (r  ,z  )  within  the  caustic  region. 

a  a  a 

In  general,  a  harmonic  component  of  a  high  frequency  wave 
pulse  satisfies 


u  (M  )  =  [A  e  10JT]M  =  g(M  ) 
a  M  a 

ci 


3u 

3z 


(M  ) 

CL 


[-iooA 


3x 
3  z 


-ioox 


e 


M 

a 


=  h(M  ) 
a 


3.7 


where  A  is  an  amplitude  function  and  only  the  term  of 
highest  order  in  frequency  has  been  kept  in  3.7.  We  assume 
that  the  direction  of  propagation  is  from  Ma  to  Mc  defined 
by  the  directional  cosines 


3t 

3r 


(M  ) 

ci 


3.8 


3  T 
3z 


(M  ) 

ci 


3t 


3  z 


a 


where  p^  and  q^  are  the  geometrical  ray  parameter  and  the 
vertical  wave  slowness  and  the  ±  signs  hold  for  downgoing 
and  upgoing  waves  respectively. 

Application  of  the  initial  conditions  3.7  to  the 
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general  solution  3.6  yields,  with  the  help  of  the  Fourier 
integral  theorem 


a±  (p) 


h 

r  Co 


a 

4tt 


00 


—CO 


g(5/Za) 


h(£,z  )  n 

ci 

icoq 


ioo  (p£±qz  ) 

ci 

e 


d£ 


Substituting  this  into  3 . 6  we  obtain 


/  r 

.  .  oj  a 

u  (r  ,  z  ,  w)  =  47  >  — 


\  h 


l 


r 


1±|1  /q  I  A(5,za)  e"1"*  d£dp 

a  J 

3.9 


where 


4>  =  ±q(z-z)  +  p (r-£)  +  t(£,z  ) 

a  a. 


3.10 


are  the  phase  functions. 

The  stationary  points  W”  =  (p  ,  £  )  of  the  phase 
functions  are  defined  by 


3  4>~ 
9P 


(PW'^W*  ° 


34>' 


(pw' V 
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This  requires 
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3t 


(P 


W'  V 


3.11 


From  3.8  it  is  apparent  that  these  conditions  are  satis¬ 
fied  if 


W1  =  (p  ,r  )  3.12 

i r  d 

i.e.,  the  geometrical  ray  parameter  and  the  initial  epi- 
central  distance  at  M  define  the  stationary  point  for  both 

d 

phase  functions.  Comparing  3.11  with  the  parametric 
equation  of  the  straight  line  ray  path  in  the  caustic 
region 


r-r 

a 


z-z 

a 


we  see  that 


|I  =  ±qr  ;  for  3.13 

a 

+ 

This  simply  means  that  $  correspond  to  upgoing  and  down¬ 
going  waves  respectively. 

For  high  frequencies,  the  double  integral  in  3 .  9  may 
be  evaluated  using  the  n-dimensional  stationary  phase  method 
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(Appendix  2).  After  simplification  using  3.7,  3.8,  3.10, 
3.12  and  3.13  the  result  is 


u  (r ,  z) 


-ij  sig  B (W) 

u  (M  )  E  (W)  e 
a 


ds 

v 


3.14 


for  both  upgoing  and  downgoing  waves .  Here 


r 


-  z-z 

1  a 

2  3 

v 


-l 


B(W)  = 


-  1 


9r 


lr  9  9 


3.15 


E(W)  = 


r  det  B(W) 


It  follows  from  3.15  that 

9r 
r3 

E(W)  = 


^2 


90 


9r 

90 


Since  cos  0  do  d $  remains  constant  along  the  ray  in  the 
caustic  region,  we  see  that  E(W)  is  in  fact  related  to  the 
ratio  of  the  ray  tube  cross-sectional  areas: 


da(r  ,z  ) 
a  a 

da (r , z) 


h 


E  (W) 


3.16 


Now  we  evaluate  the  signature  of  B  .  At  M  , |  det  B 

c 

vanishes.  This  requires  9ra/90  <  0  as  is  expected  for  a 
converging  ray  tube  radiating  outwards  from  the  source. 
Keeping  this  in  mind,  the  eigenvalues  of  the  matrix 
may  be  expressed  as 

ft •  =  h  i~  (c+d)  ±  [(c-d)2  +  4]^  f  ;  i  =  \ 

1  L  J  2 


where 


9  r  _ 

"90" 


z-z 


d  = 


2  3 

v  q 

^r 


At  the  caustic  c  =  1/d  so  that 


£L  (M  )  =  0 
1  c 


fi2(Mc)  = 


-  (c+-)  <  0 

c 


3.17 


Now 


(Mc> 


-cos  0 

0  2  3 
2v  qr 


(1±Y ) 


where 


i  = 


1 
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ranges  between  zero  and  unity.  Therefore 
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. 

7" — ~  (M  )  <  0  ;  i  =  1,2 

dr  c 


Using  this  and  3.17  it  is  clear  that  remains  negative  in 
the  vicinity  of  the  caustic  while  fi.  is  positive  for  r  <  r 

_L  C 

vanishes  at  r  =  r  and  is  negative  for  r  >  r  . 

c  c 

Therefore 


/- 


0 


r  <  r 


sig  B (W)  =  < 


3.18 


V 


-2 


r  >  r 


Using  3.16  and  3.18  in  3.14  we  have 


TT_ 


u(r,z,(o)  =  u  (M  ) 

ci 


da  (M  ) 
a 

da(r,z) 


i^H(r-rc)  -ica 


s 

ds 


v 

s 

a 


3.19 


where  H  is  the  Heavyside  function.  Therefore  we  see  there 
is  a  tt/2  phase  shift  as  a  harmonic  wave  propagates  through 
a  caustic.  The  condition  that  det  B  vanishes  at  the 
caustic  uniquely  defines  the  sign  of  all  quantities  appear¬ 
ing  in  B  so  that  3.18  always  holds  for  outgoing  waves. 
Therefore  a  tt/2  phase  shift  occurs  for  each  caustic 
encountered  along  the  ray  so  that  the  total  phase  shift 

due  ton  caustics  would  be  n  tt/2.  For  waves  propagating 
c  c 

inwards  toward  the  source,  the  same  results  can  be  obtained 
if  we  let  p  •>  -p  in  all  expressions.  It  should  be  empha¬ 
sized  that  the  sign  of  the  phase  shift,  positive  in  this 
case,  is  determined  solely  by  the  sign  convention  used  in 


. 


■ 


'■I 


the  Fourier  transform  pair.  If  the  complex  conjugate  of 
2.2  had  been  used,  the  conjugate  result,  a  -  tt/2  phase 
shift  for  each  caustic,  would  have  been  obtained. 

The  ART  result  for  the  displacement  amplitude  is 
valid  only  if  no  caustic  is  encountered  along  the  ray. 

We  now  see  that  it  is  verified  by  3.19  (r  <  r  ).  However, 
3.19  is  also  valid  for  regions  beyond  a  caustic.  Thus 
the  validity  of  the  ART  result  can  be  extended  to  these 
regions  if  a  tt/2  phase  shift  is  introduced  for  each  caustic 
encountered  along  the  ray,  i.e.,  if  we  replace  1.5  by 

TT 

my  -  1C0T 

u  (M)  =  S(0))U(M)  e  n  (M)  3.20 

where  n  ,  the  number  of  caustics,  is  given  by  3.5. 
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CHAPTER  4 

RAY  AMPLITUDE  IN  THE  VICINITY  OF  A  CAUSTIC 
A.  Introduction 

At  a  caustic,  ART  predicts  an  infinite  amplitude  for 
the  particle  displacement.  An  infinite  amplitude,  of  course, 
is  physically  meaningless  and  indicates  a  deficiency  in  the 
mathematical  formulation  of  the  theory.  To  obtain  the 
correct  dynamic  properties  of  the  wave  field  in  the  vicin¬ 
ity  of  a  caustic,  more  exact  methods  must  be  used. 

In  Chapter  2  it  was  shown  that  the  WKBJ  solution  in  the 
frequency  domain  may  be  expressed  formally  as  an  integral 
over  the  ray  parameter  p  (equation  2.20).  This  integral  may 
be  evaluated  via  different  approximation  techniques.  In 
this  chapter,  we  will  first  employ  the  second-order  saddle 
point  method  (Morse  and  Feshbach  1953,  p.  440)  to  show  that 
the  WKBJ  solution  is  equivalent  to  the  extended  ART  result 
3.20  when  the  latter  is  valid.  At  and  near  caustics  where 
the  second-order  saddle  point  method  is  inapplicable,  the 
third-order  saddle  point  method  (Chester,  Friedman  and 
Ursell  1957)  will  be  used  to  provide  the  amplitude  solu¬ 
tion.  Then,  using  a  modified  version  of  this  method,  an 
improved  solution  will  be  obtained.  Finally,  using  the 
recently  developed  equal  phase  method  (Chapman  1976,  1978), 
a  time-domain  solution  valid  over  the  entire  spatial  domain 


will  be  derived. 
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B.  Geometric  Ray  Approximation 

For  high  frequencies,  the  second-order  saddle  point 
method  may  be  applied  to  2.20  to  give 


(-1)  £S  (w)  f  (p  )  i(nj”  j[sgn(^-J)  +1]  }-iwT  (p  ,r ) 
u(M,w)  =  - ~ - r -  e  9p  pr 


where  pr  is  defined  by  the  saddle  point  condition 


9_t 

9P 


0. 


From  the  definition  of  the  phase  function  (equation  2.22) 
we  see  that  pr  must  in  fact  be  the  geometrical  ray  para¬ 
meter  of  the  ray  reaching  the  observer  point.  Therefore 
Tr  =  t  (Pr,r)  is  the  geometrical  arrival  time;  and 
the  functions  G  (equations  2.15  and  2.16)  reduce  to 


G  (M) 


(vpcosO ) 


% 

M 


for  both  P  and  SV  waves. 

Using  these  simplifications  and  3.5  for  the  number  of 
caustics  encountered  along  the  ray  we  can  write  4.1  as 

X(M,o),pr)  -io)Tr+inc  j 

U(M'U)  =  — lTM) -  e 


where 
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X(M,o>,p)  =  (-1)  S(w) 


(vp) 
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O 


(vp) 
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TT 

j=i 


(vp) 


O- 
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(vp) 


0+ 

3 


h 


R(Oj) 


and  L (M)  is  as  given  in  1.9.  This  is  called  the  geometri¬ 
cal  ray  approximation  and  is  identical  with  the  extended 
ART  result  3.20. 

C.  Airy  Approximation 

In  evaluating  the  integral  2.20  using  the  second-order 
saddle  point  method,  we  expanded  the  phase  function  t  (p,r) 
in  a  Taylor  series  about  the  saddle  point  p  and  limited 

r  2 

3  ^  T 

ourselves  to  the  quadratic  term  proportional  to  — ~  (p,r)  . 


This  approximation  is  valid  only  when 


iir 

3p2 


3d 


(P,r) 


is  not  too  small.  This  is  also  clear  from  the  fact  that 

the  geometrical  ray  approximation  4.1  tends  to  infinity  as 

.2 

d  T 


3p 


2*(p,r)  0  near  a  caustic,  which  is  physically  meaningless. 

For  regions  near  a  caustic,  we  can  expand  the  phase 


function  in  a  power  series  in  £  =  p-p  where  p  is  the 

G  C 

geometrical  ray  parameter  of  the  ray  forming  a  caustic  at 

32t 

the  epicentral  distance  r  .  Since  — j  (P  /r  )  =0,  the 

c  3  p  ^  c  c 

quadratic  term  in  the  expansion  will  be  absent.  As  a 
result 


T(p,r)  =  tc  +  +  ^''C3  +  . 


where 
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T  =  t  ( p  ,  r ) 
c  c 


T  ' 
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d  T 

9p 


<pc’r) 
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i  i  I 

c 


<pc'r) 


Using  up  to  the  cubic  term  of  this  expansion,  the  dis¬ 
placement  amplitude  2.20  may  be  approximated  by 


u  (M,  to) 


(-1)ES(m) 


,_W\  h 

K2i\r} 


,  .  /  TT 

-hot  +1  (n^-  - 
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4.2 


where 


I  = 


-ia>(x'£+  ix '  '  'K2) 
f(p)  e  c  b  c 


Bringing  the  function  f (p)  outside  the  integral  sign  at  the 

value  p  =  p  as  a  slowly  varying  function  and  introducing 
o 

a  new  variable  of  integration 


s  =  sgn  ( 1 ' ) 


4.3 


we  obtain 


I  = 
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3  f<pc> 

COT  '  '  ' 
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_ 

Y(y,s)  ds 


where 
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Y  (y ,  s) 


-i (s3/3+ys) 

=  e 


y  =  sgn  (t  '  t  '  '  ' )  ux' 
r  c  c  1  c 
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COT 


“TXT 
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4.4 


This  last  integral  is  the  Airy  function  Ai  (Abramowitz  and 
Stegun  1965)  and 


I  =  2tt 


2 


i  i  • 


COT 


1 

3  f(pc)  Ai (y ) 


Using  this  together  with  2.21  in  4.2,  we  obtain  the 
displacement  in  the  vicinity  of  a  caustic: 


u(M,co)  = 


X(M,o),pc) 

Mm} 


•  i  •  /  XT  TT « 

-lOjto+Knj-j) 


where 

L  (M)  = 


5  1 

26co^  Ai(y) 


rcosG (M)  cos0  H 

o 

^  r cosG (Ot) 

k  !  j 

TT  v(MQ)  sin0o 

i  i  cos© (Ot) 

l-i  ^  J 

Notice  that  except  when  r  =  r  ,x  is  generally  not  the 

c  c 

geometrical  arrival  time  of  a  ray;  furthermore,  the  phase 

shift  due  to  caustics  is  not  n  xr  /2  but  rather  is  the 

c 

average  of  this  value  over  the  two  geometrical  ray  branches 
forming  the  caustic. 

A  graph  of  the  Airy  function  Ai(y)  is  shown  in  Fig.  5 
It  decreases  monotonically  with  increasing  |yj  if  y  >  0, 
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Figure  5 


Graph  of  the  Airy  function  Ai . 
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and  decreases  in  an  oscillatory  manner  if  y  <  0.  The 
physical  interpretation  of  this  is  clear.  As 


T  ' 

C 
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i  i  i 
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2 


r 

c 


9p 
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y  >  0  and  y  <  0  correspond  respectively  to  the  shadow  and 
illuminated  zones  of  the  caustic.  At  each 
point  near  the  caustic  on  the  illuminated  side,  two  rays 
arrive  almost  simultaneously  and  the  oscillation  is  due  to 
the  interference  between  these  two  rays.  On  the  shadow 
side  no  rays  are  present  and  the  diffracted  energy  decays 
exponentially  with  distance  away  from  the  caustic. 

The  maximum  amplitude  does  not  occur  at  the  caustic  as 
might  be  expected  from  ART  but  is  shifted  away  from  the 
caustic  into  the  illuminated  zone.  In  general,  the  ampli¬ 
tude  peak  is  narrower  and  has  a  more  pronounced  maximum 
closer  to  the  caustic  as  the  frequency  increases. 

The  Airy  approximation  presented  here  has  been  the 
standard  method  for  the  investigation  of  the  wave  field 
near  a  caustic.  It  has  been  used  by  Brekhovskikh  (1960), 
Sachs  and  Silbiger  (1970) ,  Cerveny  and  Zahradnik  (1972) , 
Hron  and  Chapman  (1973)  ,  among  others. 


D.  Modified  Airy  Approximation 

The  airy  approximation  presupposes  that  the  function 
f (p)  (equation  2.21),  which  involves  the  product  of  reflection, 
transmission  and  surface  conversion  coefficients,  is  slowly 


•  jc 


varying  in  the  neighborhood  of  the  caustic.  We  now  pursue 
a  different  approach  which  avoids  this  limitation.  We  will 
restrict  our  consideration  to  rays  in  the  illuminated  zone 
of  the  caustic  (y  <  0) . 

Making  the  change  of  variable  4.3  but  keeping  f(s) 
inside  the  integral  we  obtain 


I  = 


1 

2 
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03T  '  '  ' 

C 

f  ( s )  Y(y,s)  ds 


4.5 


The  saddle  points  of  this  integral  are 


Sm  =  +  v/~Y 
m  JL 


m  = 


1 

2 


which  correspond  to 


p  =  p  +  sgn  (t  '  '  ' ) 
mti  ^c  c 


2x' 

c 

T  t  '  i 


1 

2 


m  = 


4.6 


Thus  it  follows  that 


r 


m  =  1,2 


Comparing  this  with  the  Taylor  expansion  of  the  geometrical 
range  about  p  =  p  ,  we  see  that  P-.  and  Pn  are  in  fact  the 
ray  parameters  of  the  two  geometrical  rays  arriving  at  the 
epicentral  distance  r.  Specifically,  we  will  see  later  that p 
and  ?2  correspond  respectively  to  the  arrival  on  the  direct 
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and  reverse  geometrical  ray  branches  which  meet  at  the 
caustic . 

We  wish  to  deform  the  path  of  integration  in  4.5  from 
the  real  s-axis  to  an  alternate  contour  in  the  complex  s- 
plane .  The  function  Y(y,s)  is  everywhere  analytic  in  the 
s-plane  and  hence  causes  no  problem.  The  function  f (s) , 
however,  involves  the  reflection  and  transmission  coeffi- 
cients  and  terms  in  q  =  (1/v^  -  p  )  2  and  therefore  requires 
more  care.  In  addition  to  the  branch  points  at  p .  =  ±  , 

J  j 

the  reflection  and  transmission  coefficients  also  have  poles 

where  the  Rayleigh  denominator  vanishes.  These  are  the 

Rayleigh  poles  which  are  always  real  and  positive 

(Archenbach  1973).  In  general,  we  see  that  all  the  poles 

of  f (p)  are  real  and  satisfy  the  inequality  |  P-i  I  ^  1/V^ ^ 

where v  is  the  maximum  velocity  encountered  at  interfaces, 
max 

The  branch  cuts  for  the  square  roots  can  be  taken  to  be  on 
the  real  p-axis  (Fig.  6).  The  saddle  points,  branch  cuts 
and  poles  in  the  complex  s-plane  are  shown  in  Fig.  7. 

Defining  the  contours  of  integration  to  be  those  as 
shown  in  Fig.  7  and  realizing  that  the  most  significant 
contributions  come  from  the  saddle  points,  we  obtain  from 


4.5. 
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If  the  endpoints  of  the  new  contours  are  kept  within  the 
shaded  zones  in  Fig.  7,  the  integrals  can  be  identified 
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Figure  6 


Poles  and  branch  cuts  on  the  p-plane  for 
the  integral  4.5. 
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Figure  7 


Poles,  branch  cuts,  saddle  points  and  contours 
on  the  s-plane  for  the  integral  4.5. 
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as  linear  combinations  of  the  Airy  functions  Ai  and  Bi 
(Appendix  3)  and 


1 

2 

3  J 

I  = 

COT  '  '  ' 

TT  < 

c 

j^f  (p-j^)  [Ai  (y)  +iBi  (y)  ] 


+  f(p2)  [Ai(y) -iBi(y) ] 


This,  together  with  2.21  and  4.2,  yields  the  displacement 
amplitude 


u  (M,  co)  =  l 
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where 
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We  denote  this  result  the  modified  Airy  approximation 


If  the  function  f (p)  is  indeed  slowly  varying  in  the 
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vicinity  of  the  caustic,  we  can  set  f(p^)  =  f  (P2 )  =  f  (Pc ) 
and  4.  7  reduces  to  the  Airy  approximation  as  is  expected. 
It  is  clear  that  the  modified  Airy  approximation  consists 
of  contributions  from  both  ray  branches.  The  linear 
combination  Ai(y)  +  iBi(y)  can  be  expressed  respectively 
as  Hankel  functions  of  one-third  order  of  the  first  and 
second  kinds  (Abramowitz  and  Stegun  (1965) )  which,  for 
| y  j  >>  1,  can  be  expanded  asymptotically  in  4 .  7  to  yield 


u  (M,co) 
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X(M,ai,pm) 

L  (M) 
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where 
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are  the  arrival  times  of  the  two  geometrical  rays. 

From  4.6  and  4.10  we  see  that  p1  always  corresponds  to 
the  ray  on  the  direct  travel  time  branch  (  3r/3p  >  0)  and 
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3  IT 

the  ray  on  the  reverse  travel  time  branch  (  <  0) . 

The  ray  on  the  reverse  branch  always  arrives  earlier  than 

the  one  on  the  direct  branch.  Thus,  we  can  conclude  that  n 

in  4 .  9  is  in  fact  the  number  of  caustics  n  (equation  3.5) 

encountered  by  the  ray.  Furthermore,  it  can  be  shown  that 

for  p  -  p  : 

^m  c 


m 


2  T  1  T  ' 

c  c 
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so  that  4.8  in  fact  are  the  geometrical  ray  amplitudes. 

Thus,  at  sufficiently  large  distances  from  the  caustic  in 
the  illuminated  zone,  the  modified  Airy  approximation  is 
equivalent  to  the  superposition  of  the  two  geometrical  ray 
amplitudes . 

Using  the  same  procedures,  the  Airy  approximation  can 
be  shown  to  consist  of  two  terms  similar  to  4.8,  except 
that  both  terms  are  evaluated  at  p  =  p  and  therefore  have 
the  same  absolute  magnitudes.  This  explains  the  complete 
destructive  interference  that  occurs  at  epicentral  distances 
where  the  Airy  function  vanishes.  More  importantly,  this 
also  shows  that  the  Airy  approximation  is  likely  to  be 
accurate  only  if  the  two  geometrical  ray  branches  have 
comparable  amplitudes  near  the  caustic.  For  wave  propagation 
in  layered  media,  this  condition  is  often  not  fulfilled. 

In  the  vicinity  of  the  PKKP  caustic,  for  instance,  the 
amplitude  of  one  ray  branch  drops  sharply  to  zero  due  to  a 
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zero  in  the  reflection  coefficient  while  the  other  branch 
remains  at  a  relatively  large  amplitude  (Richards  197  3)  .  In 
such  cases  the  modified  Airy  approximation  is  expected  to  be 
more  accurate  than  the  Airy  approximation  as  contributions 
from  the  two  ray  branches  are  given  proper  weighting. 
Numerical  results  illustrating  this  point  will  be  presented 
in  the  next  chapter . 

The  modified  Airy  approximation,  however,  is  not 

without  limitations.  Since  we  have  assumed  in  its  derivation 
3  IT 

that  -  0 ,  or  equivalently 

|r-rc|  -  0  4.11 


it  is  valid  only  in  the  neighborhood  of  a  caustic.  We  now 
determine  the  size  of  the  region  in  which  it  is  applicable. 

In  general,  the  asymptotic  expressions  for  Hankel 
functions  are  valid  if  the  value  of  the  argument  is 
comparable  with  or  greater  than  the  order  of  the  Hankel 
function.  In  our  case,  this  means  that  4.  8  follows  from 
4 .  7  provided  6  >  1  where 


3  = 


oo  I  r-r  I  2 
1  c 1 


1  9  r 

2  3  2 

9p 


2 


For  high  frequencies,  there  exists  a  region  of  epicentral 
distances  defined  by  Armj_n  <  I r  ”  rc  I  <  Armax  within 
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which  both  this  condition  and  4.11  are  satisfied 


simultaneously . 


Ar  .  is  the  value  of  r  -  r 
mm  1 


c 


at  which  3 


1  while  Ar 


max 


is  the  limit  beyond  which  4.11  becomes  a 


poor  approximation.  Within  this  region,  the  geometrical  ray 
and  modified  Airy  approximations  are  equivalent  and  both 
results  are  valid. 

For  |r  -  r  |  >  ^rmax  /  the  geometrical  ray  result  is 

certainly  applicable  as  the  caustic  is  now  even  farther 

away.  For  |r  -  r  I  <  Ar  .  ,  the  modified  Airy  approxi- 

1  c 1  mm 

mation  can  be  used  as  r  in  getting  closer  to  the  caustic 
However,  3  <  1  in  this  region  and  the  geometrical  ray 
approximation  is  invalid  as  it  is  no  longer  equivalent  to 
the  modified  Airy  approximation.  In  general,  for  high 
frequencies,  we  conclude  that  the  modified  Airy  and  geo¬ 
metrical  approximations  should  be  used  when  3  <  1  and  3 
>  1  respectively. 

E.  WKBJ  Seismogram 

We  now  consider  an  alternative  method  developed  by 
Chapman  (1976,1978)  for  the  evaluation  of  the  integral  in 
2.20.  Although  the  result  will  not  be  as  simple  for 
numerical  computation  as  those  previously  presented,  it  has 
the  advantage  of  being  valid  for  all  regions  so  that  only 
one  expression  need  be  used  at  all  times.  In  the  next 
chapter,  it  will  be  used  as  a  check  for  the  accuracy  of  the 
previous  results. 

We  invert  the  Fourier  transform  in  2. 20  to  obtain 
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u(M,t)  =  F(r,t)  * 


f(p)  6[t-T(p,r)] 


where 


F(r,t)  = 


•  /  -i  \  TT 

1  £  i(n-l)j 
e 


-rr  (2r>  h 


dS  (t)  H(t) 
dt 


^2 


and  *  denotes  a  convolution.  Convolving  this  with  a  boxcar 
function  of  width  At: 

B ( t , At)  =  [H ( t/At)  -  H (t/At-1) ] 

we  obtain  the  smoothed  time  response 


<u(M,t) > 


At 


F(r,t)  *  l 

i 


Ap. 

-AF  <f(p)>Ap. 


where 


Api  =  ipi  -  pi 


p^  and  p^  being  the  solutions  of  the  equal  phase  conditions 


t  =  x  (p ,  r ) 


and 


t  +  At  =  t  (p ,  r) 
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respectively;  and 
<f(p)>Ap.  = 


Ap± 


+ 


f (p) dp 


is  the  average  value  of  f  (p)  over  the  interval  Ap^.  This 
is  called  the  WKBJ  seismogram. 

Depending  on  whether  the  behaviour  of  the  phase  func¬ 
tion  is  essentially  quadratic  or  cubic  in  p  in  the 
neighborhood  of  interest,  we  will  see  in  the  next  chapter 
that  the  WKBJ  seismogram  is  equivalent  to  the  geometrical 
ray  and  modified  Airy  approximations.  However,  its 
numerical  evaluation  is  more  complicated  due  to  the  convo¬ 
lution  intergral. 
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CHAPTER  5 
NUMERICAL  RESULTS 

A.  Introduction 

This  chapter  presents  numerical  results  based  on  the 
extended  ART.  First,  the  various  methods  used  to  interpo¬ 
late  the  velocity  data  are  discussed.  The  corresponding 
results  for  the  geometrical  range,  travel  time,  and  the 

derivatives  and  d— are  presented  in  Appendices  4,  5, 

3p  .  ... 

and  6.  The  accuracy  of  the  Airy  and  modified  Airy  approxi¬ 
mations  is  next  compared.  Finally,  synthetic  seismograms 
for  realistic  models  are  shown. 


B.  Interpolation  of  the  Velocity  Data 

We  introduce  a  cylindrical  coordinate  system  (r,z,cf>) 
where  z,  the  depth,  is  positive  downwards.  As  before,  the 
azimuthal  angle  <j>  may  be  ignored.  Consider  a  section  of 
the  ray  specified  by  the  ray  parameter  p,  situated  at  depths 


between  z  .  and  z  .  .  If  we  denote  the  travel  time  of  this 

l  l+l 

section  of  the  ray  by  t^  and  the  corresponding  range  of 
horizontal  distance  by  r_^,  then  the  ray  integrals  are  given 


by 


t .  = 

l 


z  . 
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z  . 
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dz 

va 


r  .  = 

l 


z  . 
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pv 


dz 


where  z^  =  z^+^  if  the  ray  has  no  turning  point  and 
if  the  ray  has  a  turning  point  at  (  r|,z|  )/  and 
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a  =  (1  -  p  v  ) 


These  integrals  define  the  kinematical  properties  of  the 
ray.  To  determine  the  ray  amplitude,  we  also  must  compute 


the  derivatives  tt— 

dp 


(near  a  caustic) .  We 


denote  these  four  quantities  the  ray  characteristics. 

In  general,  the  velocity  data  consist  of  a  discrete 
set  of  velocities  at  given  depths.  The  data  must  be  inter¬ 
polated  to  provide  the  velocity  function  required  for  the 
evaluation  of  the  ray  characteristics.  The  simplest 
approach  is  to  assume  that  the  velocity  varies  linearly  be¬ 
tween  any  two  consecutive  data  points.  For  this  case, 
analytic  expressions  may  be  obtained  for  all  ray  character¬ 
istics  (Appendix  4) .  However,  piece-wise  linear  interpola¬ 
tion  introduces  false  interfaces  of  second  order  at  each 
data  point  where  the  velocity  gradient  is  discontinuous. 
This  tends  to  give  rise  to  anomalous  behaviour  of  the  ray 
amplitudes . 

To  avoid  this  difficulty,  the  velocity  data  may  be 
fitted  with  a  smoothed  spline  of  the  form 

v(z)  =  a0  +  axz  +  a2z2  +  a3z3 

(Chapman  1971)  which  ensures  the  continuity  of  the  velocity 
gradient  as  well  as  its  derivative.  This  formula  can  be 
used  insofar  as  the  data  can  be  grouped  into  sets  which 
can  be  approximated  by  piece-wise  continuous  functions. 
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With  this  interpolation,  however,  the  ray  characteristics 
must  be  evaluated  by  numerical  integration  (Appendix  5) . 

Alternatively,  we  may  apply  the  smoothed  spline 
approximation  to  interpolate  the  function  z=z (v)  instead 
of  v=v(z)  (Cerveny  1977),  i.e. 


z  (v) 


b0  +  blv  +  V2  +  b3v3 


5.1 


This  formula  is  applicable  if  the  data  can  be  grouped  into 
sets  which  can  be  approximated  by  piece-wise  monotonic 
functions.  To  compute  the  ray  characteristics  for  this 
case,  it  is  more  convenient  to  use  v  as  the  intergration 
variable : 


t . 

l 


1  dz  , 

-  -r—  dv 

va  dv 


PX  ¥  dv 
a  dv 


5.2 


where  v^=vi+1  if  the  ray  has  no  turning  point  and  v^  =  v*  if 
the  ray  has  a  turning  point  at  (r*,z|).  with  the  smoothed 
spline  interpolation  5.1,  it  is  clear  from  5.2  that  analy¬ 
tic  expressions  in  closed  form  may  be  obtained  for  the  ray 
characteristics  (Appendix  6).  This  avoids  the  necessity 
for  numerical  integration. 

Fig.  8  compares  the  primary  P  ray  seismograms  obtained 
using  the  two  different  smoothed  spline  interpolations  to 
the  velocity  function  v(z)  =  4e  Z.  The  agreement 

between  the  two  are  excellent.  However,  the  cost  in  CPU 
time  for  numerical  integration  with  v=v(z)  is  about  twice 
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Figure  8  Comparison  of  seismograms  computed  using 
numerical  integration  with  v=v(z)  and 
analytic  expressions  derived  from  z=z(v). 

The  axis  displays  the  reduced  time  where 
R  is  the  epicentral  distance  and  a 
reducing  velocity  of  8  km/sec  has  been  used. 
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that  for  the  analytic  expressions  derived  with  z  =  z (v) . 

C.  Comparisons  of  the  Airy  and  Modified  Airy  Approximations 

We  first  consider  a  simple  model  used  by  Cerveny  and 
Zahradnik  (1972)  for  the  investigation  of  the  wave  field 
near  a  caustic  (Fig.  9) .  It  consists  of  a  top  layer  with 
a  piece-wise  linear  velocity  structure  overlying  a  homo¬ 
geneous  half  space.  The  ratio  of  P  and  S  wave  velocities 
is  constant  (=/T  )  throughout  the  medium.  The  ray  diagram 
and  travel  time-distance  curves  for  primary  P  rays  are 
shown  in  Fig.  10.  Branch  1  corresponds  to  reflected  rays 
while  branches  2  and  3  represent  refracted  rays  which  form 
a  caustic  on  the  surface  at  an  epicentral  distance  of 
120  km. 

The  amplitude-distance  curves  of  the  refracted  rays 
near  the  caustic  found  by  various  methods  are  shown  in  Fig. 
11.  The  vertical  component  of  displacement  is  considered. 
The  amplitude  curves  A1  and  A 2  are  computed  using  ART.  They 
correspond  respectively  to  branches  2  and  3  of  the  travel 
time  curve  and  they  both  tend  to  infinity  at  the  caustic. 
Away  from  the  caustic,  they  are  comparable  and  remain 
essentially  constant  in  magnitude.  This  indicates  that  the 
function  f (p)  is  indeed  slowly  varying  near  the  caustic. 

A  is  the  Airy  approximation.  It  increases  from  very  small 
values  in  the  shadow  zone,  remains  finite  at  the  caustic, 
and  reaches  a  maximum  beyond  the  caustic  in  the  illuminated 
zone.  Am  is  the  modified  Airy  approximation.  It  is  almost 
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Figure  9 


Velocity  models  for  the  study  of  caustics 
the  Cerveny  and  modified  Cerveny  models. 
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VELOCITY  MODEL 
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Figure  10  Ray  diagrams  and  travel  time-distance  curves 


of  primary  P  rays  for  the  Cerveny  model. 
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Figure  11  Amplitude-distance  curves  of  refracted  P  rays 
for  the  Cerveny  model.  A^  and  A^  are  the 
geometrical  ray  amplitudes  of  the  two  ray 
branches.  The  arrows  indicate  the  direction 
of  increasing  ray  parameter.  A^  is  the 
interference  amplitude  of  these  two  branches. 
A^  is  the  Airy  approximation  and  A^  is  the 
modified  Airy  approximation. 
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identical  with  the  Airy  approximation  over  the  entire 
illuminated  zone  as  is  expected  for  a  slowly  varying  f (p) . 

A^  is  the  interference  amplitude  of  the  geometrical 
branches.  Note  that  A^,  A^  and  A  all  give  practically  the 
same  amplitudes  in  some  range  of  epicentral  distances 
beyond  the  caustic.  Immediately  near  the  caustic  ART  is 
inapplicable  and  the  Airy  and  modified  Airy  approximations 
provide  more  accurate  amplitudes. 

Next  we  consider  a  model  which  has  the  same  P  wave 

velocity  structure  as  the  previous  model.  The  S  wave 

velocity,  however,  now  differs  from  the  P  wave  velocity  by 

a  constant  amount  (v  =  v  -  3.6).  The  kinematic  proper- 

s  p 

ties  of  the  primary  P  rays,  therefore,  are  identical  with 
those  above.  The  dynamic  properties,  however,  are  now 
different . 

Fig.  12  shows  the  amplitude-distance  curves  of  the 
refracted  rays  near  the  caustic.  Note  that  the  two 
geometrical  ray  branches  A^  and  A2  no  longer  have  comparable 
amplitudes  near  the  caustic.  The  branch  A2  drops  sharply 
to  zero  just  beyond  the  caustic  due  to  a  zero  in  the  sur¬ 
face  conversion  coeffient  while  the  branch  A^  remains  at 
relatively  high  values.  Thus  f (p)  is  no  longer  slowly 
varying.  Except  at  distances  close  to  the  caustic,  the 
interference  refracted  wave  amplitude  A^  matches  closely 
the  modified  Airy  approximation  A^.  The  Airy  approximation 

A  ,  however,  yields  completely  different  results, 
a 
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. 

V  •  -  - 


Figure  12  Amplitude-distance  curves  of  refracted  P  rays 

for  the  modified  Cerveny  model.  A^  and  A ^  are 

the  geometrical  ray  amplitudes  of  the  two  ray 

branches.  The  arrows  indicate  the  direction 

of  increasing  ray  parameter.  A^  is  the 

interference  amplitude  of  these  two  branches. 

A  is  the  Airy  approximation  and  A  is  the 

m 

modified  Airy  approximation. 
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Figure  13  Comparison  of  MART  and  WKBJ  seismograms  for 
the  modified  Cerveny  model.  The  axis 
displays  the  reduced  time  where  R  is  the 
epicentral  distance  and  a  reducing  velocity 
of  6  km/sec  has  been  used. 
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To  show  that  the  modified  Airy  approximation  is  indeed 
more  accurate  than  the  Airy  approximation  in  this  case,  the 
results  are  compared  with  WKBJ  seismograms  (Fig.  13) .  On 
the  left  are  synthetic  seismograms  calculated  using  ART 
augmented  with  the  modified  Airy  approximation  near  a 
caustic  (MART  seismograms) .  The  modified  Airy  approxima¬ 
tion  is  used  for  the  MART  seismogram  at  121  km.  On  the 
right  are  the  corresponding  WKBJ  seismograms.  The  excellent 
agreement  between  the  MART  and  WKBJ  seismograms  clearly 
shows  that  the  modified  Airy  approximation  is  superior  to 
the  Airy  approximation  in  situations  where  the  product  of 
reflection,  transmission,  and  surface  conversion  coefficients 
varies  significantly  in  the  vicinity  of  the  caustic. 

D.  Synthetic  Seismograms  for  Realistic  Models 

The  FORTRAN  IV  program  developed  by  Hron  (1971)  for 
the  computation  of  synthetic  seismograms  for  layered  media 
was  extended  using  the  formulae  developed  earlier  to  include 
the  effects  of  vertical  inhomogeneity. 

Fig.  14  shows  synthetic  seismograms  for  the  vertical 
component  of  surface  displacement  for  the  Mississippi  model 
in  Fig.  15.  The  major  arrivals  are  identified  and  corre¬ 
lated  over  the  various  epicentral  distances.  Fig.  16  shows 
similar  synthetic  seismograms  for  a  modified  version  of 
the  Mississippi  model.  Here  the  inhomogeneous  layers  have 
been  replaced  by  homogeneous  layers  such  that  the  vertical 
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Figure  14  Synthetic  seismograms  for  the  Mississippi 
model.  The  axis  displays  the  reduced 
time  where  R  is  the  epicentral  distance 
and  a  reducing  velocity  of  4  km/sec  has 
been  used. 
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Figure  15  The  Mississippi  model.  The  model,  based  on 
actual  well  survey  (Narvarte  1947)  ,  is 

typical  of  certain  parts  of  Mississippi. 

The  chained  lines  show  a  modified  version 
of  the  model  where  the  inhomogeneous  layers 
are  replaced  by  homogeneous  layers  such 
that  the  vertical  travel  times  are 


preserved. 
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Figure  16  Synthetic  seismograms  for  the  homogeneous- 

layers  approximation  to  the  Mississippi  model. 
The  axis  displays  reduced  time  where  R  is 
the  epicentral  distance  and  a  reducing 
velocity  of  4  km/sec  has  been  used. 
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time  of  primary  reflections  are  preserved.  This  is  standard 
practice  in  seismic  exploration.  Thus  we  see  that  the 
general  character  of  the  seismograms  are  somewhat  similar. 
The  dynamic  properties,  however,  are  considerably  different. 
In  general,  approximation  of  an  inhomogeneous  medium  by 
homogeneous  layers  is  likely  to  be  good  only  if  a  large 
number  of  thin  layers  are  employed.  Even  then,  good 
approximation  is  attained  only  for  a  given  range  of  source 
frequencies  and  epicentral  distances  and  a  large  number  of 
rays  need  to  be  considered  (Muller  1970) .  To  study 
continuously  refracted  rays  in  media  with  continuously 
varying  properties,  it  is  much  more  efficient  to  employ 
the  extended  ART  presented  in  this  thesis. 
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Appendix  1 :  Solution  Near  a  Turning  Point 
Richards  (1974)  has  shown  that  P  and  SV  motion  in 
vertically  inhomogeneous  media  approximately  decouple  if 
the  following  potential  representation  is  used: 


u (r,z ,t) 


r 


g  (z) 


(Z)  '  v !  H(z)  °P  | 


+  ?x?xr -g(z)  $  (r,z,(|>,t)  n 

P  2(z)  b 


z 


Here  $  ,  are  the  P  and  SV  displacement  potentials  and 
the  scale  factor  g(z)  is  an  arbitrary,  smooth  function  of 
the  elastic  parameters.  At  high  frequencies,  these 
potentials  satisfy  the  wave  equations 


V  0  (r  ,  z  ,  t)  - 


V2  (z) 


3  $ (r , z,t) 

2 

at 


o 


where  v(z)  and  $  may  be  P  or  SV  wave  velocities  and 
potentials . 

Using  the  Fourier  and  Hankel  transforms  defined  in 
Chapter  2  and  assuming  azimuthal  symmetry  in  the  solution, 
we  can  reduce  this  equation  to  the  Helmholtz  equation: 


d  (p  ,  z  ,  a))  2  2*, 

- ^ -  +  a)  q  $  (p,z,co)  =  0. 

dz 


Al.l 


2 

We  assume  that  q  can  be  expanded  in  a  power  series  in  z-z* 
where  z*  is  the  turning  point,  and  that  in  the  neighbor¬ 
hood  of  z*  we  need  consider  only  the  linear  term  in  the 
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expansion.  Then 


q2  -  a(z  -  z*) 


A1.2 


where 


a 


V 


1_  dv  | 
t3  dz  *  z* 


A1.3 


This  assumption  is  invalid  only  in  the  special  cases  where 
the  turning  point  coincides  with  an  extremum  of  the  velo¬ 
city  function  and  higher  order  terms  must  be  used  (Morse 
and  Feshbach  1953,  section  9-3). 

Using  this  approximation  and  the  substitution 


1/3 


y  =  (coa) 


(z*  -  z) 


A1.4 


we  can  convert  Al.l  to  the  Stokes  equation 


dy 

Solutions  of  this  equation  are  well  known  (Budden  1961, 

Ch.  15) .  A  physically  meaningful  solution  for  the  dis¬ 
placement  potential  must  be  proportional  to  the  Airy  func¬ 
tion  Ai(y) .  If  we  consider  the  region  just  above  the 
turning  point,  i.e.  z  >  z*  or  y  <  0,  then  for  high  frequen¬ 
cies,  the  solution  has  the  form 
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where  C  is  an  arbitrary  function.  With  the  help  of  Al . 2 
to  A1.4  this  in  turn  can  be  written  as 


1 

q(5)ds  | 


z* 


This  has  the  same  form  as  the  WKBJ  solution  except  for 
the  exp(i7r/2)  phase  factor  associated  with  the  upgoing 
(totally  reflected)  wave.  Thus  the  solution  near  a  turning 
point  corresponds  to  propagation  along  the  ray  path  of  ray 
theory.  However,  analogous  to  the  problem  of  reflection 
at  interfaces,  a  phase  change  of  ir  /2  now  takes  place  as 
the  ray  is  reflected  at  the  turning  point. 
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Appendix  2 :  The  Method  of  Stationary  Phase  in  n— Dimensions 
In  this  appendix  we  derive  an  asymptotic  formula  for  a 
function  u ( w)  defined  by  the  n-dimensional  integral 


u  (a>) 


-ito$(T) 

g(T)  e  dT 

D 


A2.1 


where  T  (t^,...,  t^) ,  and  D  is  the  domain  of  integration. 
The  method  to  be  used  is  the  method  of  stationary  phase  in 
n-dimensions  (Lewis  1964) . 

A  stationary  point  of  A2 . 1  is  a  point  W  =  (w.  ,  . . . ,w  ) 

1  n 

for  which 


80 

at. 


=  o 


•  •  •  / 
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Expanding  0  about  the  stationary  point  yields 


n  n 

$  (T)  -  $  (W)  +  h\  7  (t .  -w.  )  (t .  -w  . ) 

i=i  j=i  1  1  ^  :  i: 


where 


*  =  92^(W) 

ij  “  at. at. 

i  D 


We  introduce  a  matrix  B  whose  elements  are 


B.  . 
ID 


$  .  . 
ID 


and  a  linear  coordinate  transformation  T  T'  such  that 
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the  stationary  point  is  the  origin  in  the  new  coordinates 
and  the  matrix  B  is  diagonal,  i.e.. 


$ . . (W' )  =  v . 6 . . 
ij  i  ID 


where  i=l ,  .  ..,  n  are  the  eigenvalues  of  B 

the  Kronecker  delta  function. 

In  the  new  coordinates,  then 


and  6 . . 
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A2 . 2 


where 
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a  .  =  a)  t ! 
1  1 


Using  A2 . 2  and  the  Jacobian  of  transformation 
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in  A2 . 1  we  obtain 
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n 


TT  vk 

k=l 


v,  =  det  B  (W) 


is  the  determinant  of  the  matrix  B  ;  and 


n 


I  sgn  v,  =  sig  B  (W) 


k=l 


is  the  signature  of  the  matrix  B  .  Both  these  quantities 
are  invariant  under  linear  coordinate  transformations. 
Therefore,  in  the  original  coordinates 
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Appendix  3:  Contour  Integral  Representation  of  Airy 

Functions 

Airy  functions  are  solutions  of  the  Stokes  equation 


yu 


A3.1 


We  wish  to  look  for  integral  solutions  of  the  form 


u 


•b 

a 


h(s) 


ds 


A3 . 2 


where  s  is  a  complex  variable  and  a,b  are  the  endpoints  of 
some  contour  to  be  determined.  Substitution  of  A3 . 2  into 
A3.1  and  integration  by  parts  yields 


ih(s)  e 


r 


2,  ,  ,  .  dh  :  -iys 

s  h  ( s)  -i-=—  >  e  2  ds 
ds 


0 

A3. 3 


For  the  contours  shown  in  Fig.  17,  the  first  term  of  A3 . 3 
vanishes  at  both  limits  provided  the  endpoints  of  the 
contours  are  kept  within  the  shaded  regions.  This  leads 
to  solutions  of  the  Stokes  equation  of  the  form 


u  . 
3 


Y  (y ,  s)  ds 


j  =  1*2,3 


for  some  arbitrary  constants  A^  and 

3 

-i (%r  +  ys) 

Y (y , s)  =  e 
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Figure  17  Contours  for  the  integral  A3. 3. 
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In  accordance  with  Budden  (1961)  ,  Abramowitz  and 
Stegun  (1965)  and  others,  two  independent  solutions, 
called  the  Airy  functions,  may  be  defined  as  follows: 


Ai(y) 
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2tt 


Y (y , s)  ds 


A3 .4 
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Since  Y(y,s)  is  everywhere  analytic  in  the  s-plane,  the 
contour  in  Fig.  17  can  be  deformed  to  coincide  with 
+  r2.  Then 


27TAi  (y) 
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ds  + 


Y(y,s) 


ds 


But  it  also  follows  from  A3 . 4  that 


-i27rBi(y)  = 
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Y (y , s)  ds 
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Therefore  we  obtain  the  useful  formulae 


Y(y,s)  ds  =  7T  Ai  (y)  +  iBi(y) 


Y(y,s)  ds  =  7T  Ai  (y)  -  iBi  (y)  ! 
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Appendix  4:  Ray  Characteristics  for  Piece-wise  Linear 

Interpolation 

The  velocity  function  is  locally  given  by 

v  =  v .  +  g ( z-z . ) 

1  l 

where 
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If  we  let 
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If  the  ray  segment  has  a  turning  point  at  (v^+1 , z^+1 ),  then 
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Appendix  5:  Ray  Characteristics  for  Smoothed  Spline 

Interpolation  v=v(z) 

The  velocity  function  is  locally  given  by 

0  =  b  +  b  z  +  b  z2  +  b^z3 
o  1  z  3 

where 
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z-z  . 

l 


Then,  letting  q  =  (n2-p2)^. 
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dz 


If  the  ray  segment  has  a  turning  point  at  (r.+^,z.  ), 

a  new  variable  of  integration  must  be  used  to  avoid  the 

signularity  of  z . 
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then 


r .  = 
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Appendix  6:  Ray  Characteristics  for  Smoothed  Spline 

Interpolation  z=z (V) 

The  velocity-depth  function  is  locally  given  by 


z  =  b  +  b,V  +  b  Sr2  +  b_V3 
o  1  2  3 


where  V  =  V  -  Vi  .  Then,  letting  a  =  (1-p  V  ) 
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If  the  ray  segment  has  a  turning  point  at  (r|,z|),  then 
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For  small  values  of  p,  the  inverse  square  root  in  5.2  may 
be  expanded  for  small  values  of  pv.  Keeping  only  up  to 
the  quadratic  term  of  the  expansion  the  results  are 
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where 


C0  =  bl  -  2b2vi  +  3b3Vi2 
C1  *  2(b2  -  3b3Vi> 


